################################################################################
# Created By:Pietryka
# Creation Date:  2016-08-22
# Purpose: Sensitivity analysis for CPS Turnout estimates
# Questions: mpietryka@fsu.edu
################################################################################

# PREAMBLE  ====================================



# LOAD PACKAGES --------------
library(tidyverse)



# LOAD DATA  --------------

# created in 'CPS-4A-Matching-Turnout.R'
load("Data/CPS-4A-Matching-Turnout.Rdata")


# sensitivity function
source("../Functions/sensitivity_ci.R")


#   SENSITIVITY ANALYSIS ==================================

the_row <- "Estimate"
the_col <- "prekkids"

# VALUES FOR BIAS
gamma_pos <- seq(-0.1, 0.1, 0.02)

# LABELS
label_1 <- "Assumes no bias"
label_neg <- "Assumes parents LESS likely\nto vote than non-parents"
label_pos <- "Assumes parents MORE likely\nto vote than non-parents"

# COLORS
color_1 <- "black"
color_2 <- "grey66"
color_insig <- "white"
color_annotate <- "grey40"

# 2010  --------------

estimates <- tibble(
  year = c(2010L, 2014L),
  est = list(turnout_est_10, turnout_est_14),
  mod = map(est, "att.model"),
  beta = map_dbl(mod, ~ .x["Estimate", "prekkids"]),
  se   = map_dbl(mod, ~ .x["Std. Error", "prekkids"]),
  lb   = beta - 1.96 * se,
  ub   = beta + 1.96 * se)   %>%
  select(-mod) %>%
  mutate(sens_2 = map2(beta, se, ~
                         s_ci(.x, .y,
                                 gamma = gamma_pos,
                                 delta = 0.2,
                                 mod_type = "lm")))  %>%
  mutate(sens_4 = map2(beta, se, ~
                         s_ci(.x, .y,
                                 gamma = gamma_pos,
                                 delta = 0.4,
                                 mod_type = "lm")))  %>%
  mutate(sens_6 = map2(beta, se, ~
                         s_ci(.x, .y,
                                 gamma = gamma_pos,
                                 delta = 0.6,
                                 mod_type = "lm")))  %>%
  gather(delta, sens, starts_with("sens_"))  %>%
  mutate(delta = as.numeric(gsub("sens_", "", delta))/10)  %>%
  mutate(delta = paste("delta==", delta))  %>%
  mutate(gamma = map(sens, "gamma"))  %>%
  mutate(beta_crt = map(sens, "corrected_beta"))  %>%
  mutate(lb_crt = map(sens, "corrected_lb"))  %>%
  mutate(ub_crt = map(sens, "corrected_ub"))   %>%
  select(-sens, -est)  %>%
  unnest()   %>%
  mutate(dot_type = case_when(
    beta == beta_crt ~ "original",
    sign(lb_crt) != sign(ub_crt) ~ "statistically insignificant",
    beta != beta_crt ~ "corrected"
  ))  %>%
  mutate(line_type = case_when(
    beta == beta_crt ~ "original",
    beta != beta_crt ~ "corrected"
  ))  %>%
  mutate(similarity = beta - beta_crt)  %>%
  mutate(the_label =  recode(round(gamma, 1),
                             `0.9` = label_neg,
                             `1.0` = label_1,
                             `1.1` = label_pos,
                             .default = "")) %>%
  mutate(the_color = ifelse(gamma == 0, color_1, color_2))  %>%
  mutate(sig_color = ifelse(sign(lb_crt) == sign(ub_crt),
                            the_color,
                            color_insig))





# Figure A1 =============================================



estimates_10 <- estimates  %>%
  filter(year == 2010 & delta == "delta== 0.4")
estimates_14 <- estimates  %>%
  filter(year == 2014 & delta == "delta== 0.4")

label_data <- tibble(
  the_label = c(label_neg, label_1, label_pos),
  the_x     = c(-.02, 0, .02),
  the_y     = c(-0.133, -0.08, -0.133),
  the_col   = c(color_annotate, color_1, color_annotate)
)

arrow_y <- -0.12
arrow_start_x <- 0.02



or_plot <- function(data, year = "", add_labs = FALSE,  the_cap = "") {
  the_plot <- ggplot(data, aes(x = gamma,
                               y = beta_crt,
                               label = the_label)) +
    geom_hline(yintercept = 0, linetype = 2, size = 1.5) +
    geom_pointrange(aes(ymin = lb_crt, ymax = ub_crt),
                    shape = 21,
                    fatten = 4,
                    linewidth = 1.1,
                    color = data$the_color,
                    fill =   data$sig_color) +
    #scale_y_continuous(breaks = seq(.6, 1.3, .2)) +
    expand_limits(y = c(-.14, 0.05)) +
    xlab(    expression(gamma * " (Difference in means)")) +
    ylab("") +
    labs(title = year,
         caption = the_cap) +
    #coord_flip() +
    theme_minimal() +
    theme(plot.caption = element_text(colour = color_annotate)) +
    theme(legend.position = "none") +
    theme(text = element_text(size = 16))
  if(add_labs == TRUE){
    the_plot <- the_plot +
      geom_text(data = label_data,
                aes(x = the_x, y = the_y, label = the_label),
                hjust = "outward",
                size = 3.5,
                color = label_data$the_col,
                fontface = "bold") +
      geom_segment(x = arrow_start_x,
                   xend = gamma_pos[11],
                   y = arrow_y,
                   yend = arrow_y,
                   arrow = arrow(),
                   color = color_annotate,
                   size = 2) +
      geom_segment(x = arrow_start_x*-1,
                   xend = gamma_pos[1],
                   y = arrow_y,
                   yend = arrow_y,
                   arrow = arrow(),
                   color = color_annotate,
                   size = 2) +
      ylab("Effect of Parenthood on Self-Reported Turnout\n(Difference in means)")
  }
  the_plot
}


graph_10 <- or_plot(estimates_10, year = 2010, add_labs = TRUE, the_cap = "")
graph_14 <- or_plot(estimates_14,
                    year = 2014,
                    the_cap = NULL)



# Figure A2 =============================================




the_plot <- ggplot(estimates, aes(x = gamma,
                      y = beta_crt,
                      color = line_type,
                      fill = dot_type)) +
  facet_grid(year ~ delta, scales = "free_y",
             labeller = labeller(.cols = label_parsed)) +
  geom_hline(yintercept = 0, linetype = 2, size = 1) +
  geom_pointrange(
    aes(ymin = lb_crt, ymax = ub_crt),
    shape = 21,
    alpha = 1,
    fatten = 3,
    linewidth = 1
  ) +
  xlab(expression(gamma * ", effect of the omitted variable on self-reported turnout (difference in means)")) +
  ylab("Effect of parenthood on self-reported turnout\n(difference in means)") +
  expand_limits(x = c(-.11, .11)) +
  theme_minimal(base_size = 12) +
  theme(plot.caption = element_text(colour = color_annotate)) +
  theme(legend.position = "none") +
  theme(panel.grid.minor = element_blank()) +
  scale_color_manual(values = c("grey66", "black")) +
  scale_fill_manual(values = c("grey66", "black", "white"))




# Save =======================================

graphics.off()
windows(12, 6)
ggpubr::ggarrange(graph_10, graph_14, ncol = 2, align = "h")
ggsave("Figures/CPS-Turnout-Sensitivity-Main.tiff", compression = "lzw")


graphics.off()
windows(6.5, 8)
the_plot
ggsave("Figures/CPS-Turnout-Sensitivity-Appendix.tiff", compression = "lzw")


